Appendix - Code
#setup and cleaning data
library("tidyverse")
library("scales")
library("MASS")
library("viridis")
library("DT")
library("rcompanion")
library("ggpubr")
library("car")
library("e1071")
library("gridExtra")
library("reshape2")
library("lubridate")
setwd("/Users/edenaxelrad/Desktop/Green Roofs (Thesis)/Summer Run")
treat_avg <- read.csv("Interior Avgs.csv", stringsAsFactors = FALSE)
treat_avg$timestamp <- mdy_hm(treat_avg$Date.Time)
treat_avg$hour <- hour(treat_avg$timestamp)
treat_avg$day <- date(treat_avg$timestamp)
#Background information for treatment averaged data
head <- head(treat_avg)
knitr::kable(head, caption = "Data Headers and Sample Records")
rubro <- cbind.data.frame(treat_avg$timestamp, treat_avg$Sedum.rubro, "Sedum rubrotinctum")
colnames(rubro) <- c("Date", "Temperature", "Treatment")
acre <- cbind.data.frame(treat_avg$timestamp, treat_avg$Sedum.acre, "Sedum acre")
colnames(acre) <- c("Date", "Temperature", "Treatment")
soil <- cbind.data.frame(treat_avg$timestamp, treat_avg$bare.substrate, "Bare substrate")
colnames(soil) <- c("Date", "Temperature", "Treatment")
treat_temps <- rbind.data.frame(rubro, acre, soil)
ggplot(treat_temps, aes(Date, Temperature, color = Treatment)) +
geom_path() + facet_grid(Treatment ~ .) +
theme(legend.position = "none") +
ggtitle("Temperature by Treatment")
# Comparison of Treatment Avg at Samples
rubro <- cbind(treat_avg$rubro.vs.substrate, "rubro vs substrate")
colnames(rubro) <- c("difference", "comparison")
acre <- cbind(treat_avg$acre.vs.substrate, "acre vs substrate")
colnames(acre) <- c("difference", "comparison")
plot(treat_avg$Sample_Number, treat_avg$rubro.vs.substrate, type = "l",
main = "(a) sedum rubrotinctum - bare substrate",
xlab = "sample number",
ylab = "temperature difference (C)")
abline(h=0, col="blue")
abline(h=mean(treat_avg$rubro.vs.substrate), col="red", lty = 2)
plot(treat_avg$Sample_Number, treat_avg$acre.vs.substrate, type = "l",
main = "(b) sedum acre - bare substrate",
xlab = "sample number",
ylab = "temperature difference (C)")
abline(h=0, col="blue")
abline(h=mean(treat_avg$acre.vs.substrate), col="red", lty = 2)
long_data <- rbind.data.frame(rubro, acre)
long_data$difference <- as.numeric(as.character(long_data$difference))
long_data$abs_diff <- lapply(long_data$difference, abs)
summ_stats <-
long_data %>%
group_by(comparison) %>%
summarize(mean(difference), max(difference), sd(difference))
knitr::kable(summ_stats, digits = 4, caption = "Table: Summary Statistics")
# Temperature Differences by Hour
neg_counts <- c("negative/colder", sum(treat_avg$rubro.vs.substrate<0), sum(treat_avg$acre.vs.substrate<0))
pos_counts <- c("positive/hotter", sum(treat_avg$rubro.vs.substrate>0), sum(treat_avg$acre.vs.substrate>0))
zero_counts <- c("zero/same", sum(treat_avg$rubro.vs.substrate==0), sum(treat_avg$acre.vs.substrate==0))
diff_counts <- rbind.data.frame(neg_counts, pos_counts, zero_counts)
colnames(diff_counts) <- c("counts", "rubro - substrate", "acre - substrate")
diff_counts$`rubro - substrate` <- as.numeric(as.character(diff_counts$`rubro - substrate`))
diff_counts$`acre - substrate` <- as.numeric(as.character(diff_counts$`acre - substrate`))
diff_counts <-
diff_counts %>%
mutate("Percent (%) rubro" = signif(100*(`rubro - substrate`/1589), digits=5)) %>%
mutate("Percent (%) acre" = signif(100*(`acre - substrate`/1589), digits=4))
diff_countsR <- diff_counts[,c(1,2,4)]
diff_countsA <- diff_counts[,c(1,3,5)]
get_density <- function(x, y, ...) {
dens <- MASS::kde2d(x, y, ...)
ix <- findInterval(x, dens$x)
iy <- findInterval(y, dens$y)
ii <- cbind(ix, iy)
return(dens$z[ii])
}
treat_avg$hour <- as.numeric(treat_avg$hour)
treat_avg$densityR <- get_density(treat_avg$hour, treat_avg$rubro.vs.substrate, n=40)
treat_avg$densityA <- get_density(treat_avg$hour, treat_avg$acre.vs.substrate, n=40)
ggplot(treat_avg) +
geom_point(aes(treat_avg$hour, treat_avg$rubro.vs.substrate, size = treat_avg$densityR,
color = abs(treat_avg$rubro.vs.substrate))) +
geom_hline(yintercept=0, color = "white", alpha = 0.9) +
geom_smooth(aes(treat_avg$hour, treat_avg$rubro.vs.substrate), color = "darkblue", se = FALSE) +
scale_color_continuous(low = "blue", high = "red") +
guides(color = FALSE, size = FALSE) +
ylim(-3.5,7.5) + ggtitle("(a) Sedum rubro - Bare Substrate") +
ylab("Temperature Difference") +
xlab("Hour")
knitr::kable(diff_countsR, digits = 4, caption = "Table: Difference Score Counts for Sedum rubrotinctum")
ggplot(treat_avg) +
geom_point(aes(treat_avg$hour, treat_avg$acre.vs.substrate, size = treat_avg$densityA,
color = abs(treat_avg$acre.vs.substrate))) +
geom_hline(yintercept=0, color = "white", alpha = 0.9) +
geom_smooth(aes(treat_avg$hour, treat_avg$acre.vs.substrate), color = "darkblue", se = FALSE) +
scale_color_continuous(low = "blue", high = "red") +
guides(color = FALSE, size = FALSE) +
ylim(-3.5,7.5) + ggtitle("(b) Sedum acre - Bare Substrate") +
ylab("Temperature Difference") +
xlab("Hour")
knitr::kable(diff_countsA, digits = 4, caption = "Table: Difference Score Counts for Sedum acre")
# Comparison of Temperatures by Hour
hrly_avg <-
treat_avg %>%
group_by(hour) %>%
summarize(mean(Sedum.rubro), mean(Sedum.acre), mean(bare.substrate))
hrly_avg <- hrly_avg[1:24,]
hrly_avg$rubro_diff <- hrly_avg$`mean(Sedum.rubro)`- hrly_avg$`mean(bare.substrate)`
hrly_avg$acre_diff <- hrly_avg$`mean(Sedum.acre)`- hrly_avg$`mean(bare.substrate)`
hrly_avg[,2:6] <- round(hrly_avg[,2:6], 2)
plot(hrly_avg$hour, hrly_avg$`mean(Sedum.acre)`, type = "l", col = "blue",
main = "Hourly Temperature Averages",
xlab = "hour",
ylim = c(15,40),
ylab = "avg temperature (C)")
lines(hrly_avg$hour, hrly_avg$`mean(bare.substrate)`, col="black")
lines(hrly_avg$hour, hrly_avg$`mean(Sedum.rubro)`, col="red")
legend("topleft", legend = c("substrate", "rubro", "acre"), col = c("black", "red", "blue"), lty=1:1, box.lty=0)
only_8am <-
treat_avg %>%
filter(hour == "8" | hour == "08") %>%
group_by(day) %>%
summarize(mean(bare.substrate), mean(Sedum.rubro), mean(Sedum.acre))
plot(only_8am$day, only_8am$`mean(bare.substrate)`, type = "l", main = "Daily Temperatures at 8am", xlab = "", ylab = "Temperature (C)", ylim = c(10,25))
lines(only_8am$day, only_8am$`mean(Sedum.rubro)`, type = "l", col = "red")
lines(only_8am$day, only_8am$`mean(Sedum.acre)`, type = "l", col = "blue")
legend("topleft", legend = c("substrate", "rubro", "acre"), col = c("black", "red", "blue"), lty=1:1, box.lty=0)
only_daytime <-
treat_avg %>%
filter(hour > "10") %>%
filter(hour < "17") %>%
group_by(day) %>%
summarize(mean(bare.substrate), mean(Sedum.rubro), mean(Sedum.acre))
plot(only_daytime$day, only_daytime$`mean(bare.substrate)`, type = "l", main = "Daily Temperatures from 11am to 4pm", xlab = "", ylab = "Temperature (C)", ylim = c(25, 42))
lines(only_daytime$day, only_daytime$`mean(Sedum.rubro)`, type = "l", col = "red")
lines(only_daytime$day, only_daytime$`mean(Sedum.acre)`, type = "l", col = "blue")
legend("topleft", legend = c("substrate", "rubro", "acre"), col = c("black", "red", "blue"), lty=1:1, box.lty=0)
datatable(hrly_avg, class = 'cell-border stripe', options = list(pageLength = 24))
# Diurnal Temperature Ranges (DTR)
daily_range <-
treat_avg %>%
group_by(day) %>%
summarize(max(Sedum.rubro), min(Sedum.rubro), mean(Sedum.rubro), max(Sedum.acre), min(Sedum.acre), mean(Sedum.acre), max(bare.substrate), min(bare.substrate), mean(bare.substrate))
daily_range <-
daily_range %>%
mutate(rubro_range = `max(Sedum.rubro)` - `min(Sedum.rubro)`) %>%
mutate(acre_range = `max(Sedum.acre)` - `min(Sedum.acre)`) %>%
mutate(soil_range = `max(bare.substrate)` - `min(bare.substrate)`)
plot(daily_range$day, daily_range$rubro_range, type = "l", col = "red",
main = "Daily Temperature Ranges (Interior Sensors)",
xlab = "",
ylab = "Temperature Range (C)",
ylim = c(10,25))
lines(daily_range$day, daily_range$acre_range, col="blue", lty = 1)
lines(daily_range$day, daily_range$soil_range, col="black", lty = 1)
legend("top", legend = c("substrate", "rubro", "acre"), col = c("black", "red", "blue"), lty=1:1, box.lty=0)
DTR_summ <-
daily_range %>%
summarize(mean(rubro_range), mean(acre_range), mean(soil_range))
knitr::kable(DTR_summ, digits = 4, caption = "Table: Summary of Mean DTR for Each Treatment")
Count_types <- c("Smaller DTR", "Larger DTR", "Same DTR")
Rdtr_counts <- rbind(sum((daily_range$rubro_range - daily_range$soil_range) < 0), sum((daily_range$rubro_range - daily_range$soil_range) > 0), sum((daily_range$rubro_range - daily_range$soil_range) == 0))
Adtr_counts <- rbind(sum((daily_range$acre_range - daily_range$soil_range) < 0), sum((daily_range$acre_range - daily_range$soil_range) > 0), sum((daily_range$acre_range - daily_range$soil_range) == 0))
DTR_counts <- cbind.data.frame(Count_types, Rdtr_counts, Adtr_counts)
# Daily Max Temperatures
plot(daily_range$day, daily_range$`max(Sedum.rubro)`, type = "l", col = "red",
main = "Daily Max Temperature (Interior Sensors)",
xlab = "",
ylab = "Temperature (C)",
ylim = c(30,42))
lines(daily_range$day, daily_range$`max(Sedum.acre)`, col="blue", lty = 1)
lines(daily_range$day, daily_range$`max(bare.substrate)`, col="black", lty = 1)
legend("top", legend = c("substrate", "rubro", "acre"), col = c("black", "red", "blue"), lty=1:1, box.lty=0)
Category <- c("Lower Max (vs Soil)", "Higher Max (vs Soil)", "Same Max (vs Soil)", "Average Max Temperature", "Maximum Max Temperature")
rubro_max <- rbind(sum((daily_range$`max(Sedum.rubro)` - daily_range$`max(bare.substrate)`) < 0),
sum((daily_range$`max(Sedum.rubro)` - daily_range$`max(bare.substrate)`) > 0),
sum((daily_range$`max(Sedum.rubro)` - daily_range$`max(bare.substrate)`) == 0),
mean(daily_range$`max(Sedum.rubro)`),
max(daily_range$`max(Sedum.rubro)`))
acre_max <- rbind(sum((daily_range$`max(Sedum.acre)` - daily_range$`max(bare.substrate)`) < 0),
sum((daily_range$`max(Sedum.acre)` - daily_range$`max(bare.substrate)`) > 0),
sum((daily_range$`max(Sedum.acre)` - daily_range$`max(bare.substrate)`) == 0),
mean(daily_range$`max(Sedum.acre)`),
max(daily_range$`max(Sedum.acre)`))
Max_counts <- cbind.data.frame(Category, rubro_max, acre_max)
knitr::kable(Max_counts, digits = 3, caption = "Table: Summary of Maximum Temperature by Treatment")
# Daily Min Temperatures
plot(daily_range$day, daily_range$`min(Sedum.rubro)`, type = "l", col = "red",
main = "Daily Min Temperature (Interior Sensors)",
xlab = "",
ylab = "Temperature (C)",
ylim = c(12,25))
lines(daily_range$day, daily_range$`min(Sedum.acre)`, col="blue", lty = 1)
lines(daily_range$day, daily_range$`min(bare.substrate)`, col="black", lty = 1)
legend("top", legend = c("substrate", "rubro", "acre"), col = c("black", "red", "blue"), lty=1:1, box.lty=0)
Category <- c("Lower Min (vs Soil)", "Higher Min (vs Soil)", "Same Min (vs Soil)", "Average Min Temperature", "Minimum Min Temperature")
rubro_max <- rbind(sum((daily_range$`min(Sedum.rubro)` - daily_range$`min(bare.substrate)`) < 0),
sum((daily_range$`min(Sedum.rubro)` - daily_range$`min(bare.substrate)`) > 0),
sum((daily_range$`min(Sedum.rubro)` - daily_range$`min(bare.substrate)`) == 0),
mean(daily_range$`min(Sedum.rubro)`),
min(daily_range$`min(Sedum.rubro)`))
acre_max <- rbind(sum((daily_range$`min(Sedum.acre)` - daily_range$`min(bare.substrate)`) < 0),
sum((daily_range$`min(Sedum.acre)` - daily_range$`min(bare.substrate)`) > 0),
sum((daily_range$`min(Sedum.acre)` - daily_range$`min(bare.substrate)`) == 0),
mean(daily_range$`min(Sedum.acre)`),
min(daily_range$`min(Sedum.acre)`))
Max_counts <- cbind.data.frame(Category, rubro_max, acre_max)
knitr::kable(Max_counts, digits = 3, caption = "Table: Summary of Minimum Temperature by Treatment")
t_test <- c("rubro DTR", "acre DTR", "rubro max", "acre max", "rubro min", "acre min")
p_value <- c(t.test(daily_range$rubro_range, daily_range$soil_range, paired = TRUE)$p.value,
t.test(daily_range$acre_range, daily_range$soil_range, paired = TRUE)$p.value,
t.test(daily_range$`max(Sedum.rubro)`, daily_range$`max(bare.substrate)`, paired = TRUE)$p.value,
t.test(daily_range$`max(Sedum.acre)`, daily_range$`max(bare.substrate)`, paired = TRUE)$p.value,
t.test(daily_range$`min(Sedum.rubro)`, daily_range$`min(bare.substrate)`, paired = TRUE)$p.value,
t.test(daily_range$`min(Sedum.acre)`, daily_range$`min(bare.substrate)`, paired = TRUE)$p.value)
CI_Lower <- c(t.test(daily_range$rubro_range, daily_range$soil_range, paired = TRUE)$conf.int[1],
t.test(daily_range$acre_range, daily_range$soil_range, paired = TRUE)$conf.int[1],
t.test(daily_range$`max(Sedum.rubro)`, daily_range$`max(bare.substrate)`, paired=TRUE)$conf.int[1],
t.test(daily_range$`max(Sedum.acre)`, daily_range$`max(bare.substrate)`, paired=TRUE)$conf.int[1],
t.test(daily_range$`min(Sedum.rubro)`, daily_range$`min(bare.substrate)`,paired=TRUE)$conf.int[1],
t.test(daily_range$`min(Sedum.acre)`, daily_range$`min(bare.substrate)`, paired=TRUE)$conf.int[1])
CI_Upper <- c(t.test(daily_range$rubro_range, daily_range$soil_range, paired = TRUE)$conf.int[2],
t.test(daily_range$acre_range, daily_range$soil_range, paired = TRUE)$conf.int[2],
t.test(daily_range$`max(Sedum.rubro)`, daily_range$`max(bare.substrate)`, paired=TRUE)$conf.int[2],
t.test(daily_range$`max(Sedum.acre)`, daily_range$`max(bare.substrate)`, paired=TRUE)$conf.int[2],
t.test(daily_range$`min(Sedum.rubro)`, daily_range$`min(bare.substrate)`,paired=TRUE)$conf.int[2],
t.test(daily_range$`min(Sedum.acre)`, daily_range$`min(bare.substrate)`, paired=TRUE)$conf.int[2])
ttest_pvals <- cbind.data.frame(t_test, p_value, CI_Lower, CI_Upper)
knitr::kable(ttest_pvals, digits = 4, caption = "Table: Summary of Paired T-tests for DTR, Min, and Max Temperatures (alpha = 0.05)")
# Linear Models for Treatment vs Control
ggplot(treat_avg, aes(x=bare.substrate, y=Sedum.rubro)) +
geom_point()+
geom_smooth(method = "lm")+
ylim(10,43) + xlim(10,43) +
ylab("sedum rubro")+
xlab("bare substrate")+
geom_abline(b=1, color="red", linetype=2)+
ggtitle("(a) Linear Model of soil ~ Sedum rubrotinctum")
linearMod1 <- lm(bare.substrate ~ Sedum.rubro , data=treat_avg) %>%
summary()%>%
broom::tidy()
knitr::kable(linearMod1, digits = 4, caption = "Linear Model Summary Statistics: soil ~ Sedum rubrotinctum")
ggplot(treat_avg, aes(x=bare.substrate, y=Sedum.acre)) +
geom_point()+
geom_smooth(method = "lm")+
ylim(10,43) + xlim(10,43) +
ylab("sedum acre")+
xlab("bare substrate")+
geom_abline(b=1, color="red", linetype=2)+
ggtitle("(b) Linear Model of soil ~ Sedum acre")
linearMod2 <- lm(bare.substrate ~ Sedum.acre , data=treat_avg) %>%
summary()%>%
broom::tidy()
knitr::kable(linearMod2, digits = 4, caption = "Linear Model Summary Statistics: soil ~ Sedum acre")
ggplot(treat_avg, aes(x=bare.substrate, y=rubro.vs.substrate)) +
geom_point()+
geom_smooth()+
ylab("sedum rubro. - bare substrate")+
xlab("bare substrate")+
ggtitle("(a) Loess regression of soil ~ (Sedum rubrotinctum - soil)")
linearMod3 <- lm(bare.substrate ~ rubro.vs.substrate, data=treat_avg) %>%
summary() %>%
broom::tidy()
knitr::kable(linearMod3, digits = 4, caption = "Linear Model Summary Statistics")
ggplot(treat_avg, aes(x=bare.substrate, y=acre.vs.substrate)) +
geom_point()+
geom_smooth()+
ylab("sedum acre - bare substrate")+
xlab("bare substrate")+
ggtitle("(b) Loess regression of soil ~ (Sedum acre - soil)")
linearMod4 <- lm(bare.substrate ~ acre.vs.substrate, data=treat_avg) %>%
summary() %>%
broom::tidy()
knitr::kable(linearMod4, digits = 4, caption = "Linear Model Summary Statistics")
#looking at loess for daytime only - 10am to 7pm
daytime <- treat_avg %>%
filter(hour > "10") %>%
filter(hour < "19")
ggplot(daytime, aes(x=bare.substrate, y=acre.vs.substrate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_point()+
geom_smooth(method = "lm")+
xlim(20,45) + ylim(-5,5) +
ylab("sedum acre - bare substrate") +
xlab("bare substrate") +
ggtitle("LM of soil ~ (Sedum acre - soil) during the day time (10am to 7pm)")
# Appendix - Code